

#------------------------------------------Housekeeping-------------------------------------------------------------#
#
# This is the primary replication file for "How windfall income increases spending at poker machines". 
# 
# This code replicates all analyses with one exception: 
# "Figure 1: Average real monthly expenditures per EGM" which is created by LGAmapsV3.R 
#
# created 21 Sept. 2012 by Hielke Buddelmeyer
# last update 17 Jan. 2013 by Kyle Peyton
# 
# Questions, comments, complaints? contact Kyle Peyton: kyle.peyton@unimelb.edu.au 
#-------------------------------------------------------------------------------------------------------------------#

rm(list=ls())

#Packages needed to reproduce this analysis. You will need to download these if not already installed. 
library('ggplot2')
library('xtable')
library('foreign')
library('Zelig')
library('scales')
library('pcse')


#------------------------------------------Data management-------------------------------------------#
#                                                                                                    #
#                                                                                                    #
#        Create the estimated monthly anomoly (EMA) dummies: 96 of these, one for each period        #          
#        (=specific month-year combinations)                                                         #
#                                                                                                    # 
#----------------------------------------------------------------------------------------------------#

#import pokies data; path name depends on your system. 
pokies <- read.dta('C:/Users/peytonk/Dropbox/Pokies/replication/pokies01.dta') #PC

#create factor variables necessary for fixed effects specification. First for LGA, next for time.
pokies$LGA <- factor(pokies$LGA_short)
table(pokies$LGA)

pokies$time <- factor(pokies$Period)
table(pokies$time)

#Generate 96 EMA dummies, one for each period 
pokies$PeriodCopy <- pokies$Period 
for(i in 1:96){ 
  pokies[[paste("EMA", i, sep = "_")]] <- (pokies$PeriodCopy==i)*1
  }

#----------------------------------------Descriptive statistics------------------------------------------#
#
#  Time series plots to appear in the Appendix
# 
# Note: output directories for pdfs depends on your system! 
#--------------------------------------------------------------------------------------------------------#

# Trim off excess margin space (bottom, left, top, right)
par(mar=c(2.5, 2.5, 1.5, 0.2))

#Decomposition of additive time series for Victoria (Figure 4 in paper)
victotal <- as.matrix(subset(pokies, LGA_short == "Vic_Total",select = c(rMonExp_EGM)))
tsvic <- ts(victotal/100, frequency = 12, start=c(2004,7))

pdf(file='C:/Users/peytonk/Dropbox/Pokies/Replication/vicdecomp.pdf', height = 8.5, width = 11)
tsviccomponents <- decompose(tsvic, type = "additive")
plot(tsviccomponents, type = "b", xlab = "", xlim = c(2004,2012))
dev.off()

#Decomposition of additive time series for Surf Coast (Figure 5 in paper)
surf <- as.matrix(subset(pokies, LGA_short == 'Surf_Coast', select = c(rMonExp_EGM)))
tssurf <- ts(surf/100, frequency = 12, start=c(2004,7))
plot.ts(tssurf)

pdf(file='C:/Users/peytonk/Dropbox/Pokies/Replication/surfedecomp.pdf', height = 8.5, width = 11)
tssurfcomponents <- decompose(tssurf)
plot(tssurfcomponents, type = "b", xlab = "", xlim = c(2004,2012))
dev.off()


#--------------------------------------------SIMPLE MODELS-------------------------------------------------#
#
# Used for descriptive purposes in paper. No tables presented from these results. 
#
#-----------------------------------------------------------------------------------------------------------#

#Linear regression without LGA fixed effects
z.out1 <- zelig(rMonExp_EGM~Month2+Month3+Month4+Month5+Month6+Month7+Month8+Month9+Month10+Month11
                +Month12+Yrs_minus_2004, model="ls", data=pokies)


#Linear regression with LGA fixed effects
z.out2 <- zelig(rMonExp_EGM~Month2+Month3+Month4+Month5+Month6+Month7+Month8+Month9+Month10+Month11
               +Month12+Yrs_minus_2004 + as.factor(LGA), model="ls", data=pokies)

summary(z.out1)
summary(z.out2)

#------------------------------EMA regression----------------------------------------------# 
# This corresponds to equation 4 in the paper. The results are used to produce Figure 2 
# and Table 4. 
#------------------------------------------------------------------------------------------#

#Increase memory to 3GB; this will not work on a Linux machine 
memory.limit(size=300000)

#We have 63 lga variables because LGA_short54 is "Vic_Total", which we exclude 
ema_variables <- paste('EMA', 1:96, sep="_")

#original regression formula: real monthly expenditures per egm on monthly dummies using 'Jan' as 
#the base + a linear time trend + LGA fixed effects 
ema_regression <- rMonExp_EGM~ Month2+Month3+Month4+Month5+Month6+Month7+Month8+Month9+
  Month10+Month11+Month12+Yrs_minus_2004 + as.factor(LGA) - 1

#Create list of 96 copies of the original regression
ema_regression_list <- lapply(ema_variables, function(x, orig = ema_regression) { 
    new <- reformulate(c(x,'.'))
    update(orig, new)})

#apply the linear model to the list of 96 regression fomulas using the pokie data 
models <- lapply(ema_regression_list, lm, data = pokies)

#name the models 'reg.out1' to 'reg.out96' 
names(models) <-paste('reg.out', 1:96, sep="")

#create summarries for all the models
summaries <- lapply(models, summary) 

#Extract the coefficients from the models 
coefficients <- lapply(models, coef)

#Create table with coefficient estimates and standard errors for all 96 models
coef_tables <- lapply(summaries, '[[', 'coefficients')

#Extract 96 coefficients for estimated monthly anomaly (EMA) 
ema.beta <- matrix(nrow=96,ncol=1) 

for (i in 1:96){
  ema.beta[i,] <- c(models[[paste('reg.out',i,sep='')]]$coefficients[1])
}

#Extract 96 classical standard errors for the coefficients 
ema.se <- matrix(nrow=96,ncol=1) 

for (i in 1:96){
  ema.se[i,] = c(sqrt(diag(vcov(models[[paste('reg.out',i,sep='')]])))[1])
}

#Extract 96 panel corrected standard errors using pcse (this takes some time). 
ema.pcse <- matrix(nrow=96,ncol=1) 

for (i in 1:96){
  ema.pcse[i,] = c(pcse(models[[paste('reg.out',i,sep='')]], groupN = pokies$LGA, groupT = pokies$Period)$pcse[1])
}

#-------------------------------------------------StarGazing-------------------------------------------------------# 
#
# Create Table 4 from the paper: estimates of monthly anomalies from equation 4   
#
# Note: we have used panel corrected standard errors because we think it's the most honest way to present the results. 
# You are free to use classical errors, robust, etc. Whatever tickles your methodological fancy. Note that the plain
# vanilla OLS estimates are only technically valid when the off-diagonals of the variance matrix from the regression
# are equal to zero (independence) and 'robust' or White's standard errors only address heteroscedasticity. 
#------------------------------------------------------------------------------------------------------------------#

#create data frame of betas and pcses for export to latex 
stargaze = as.data.frame(cbind(ema.beta, ema.pcse))

starlab = c("Jul","Aug","Sep","Oct","Nov","Dec", "Jan","Feb","Mar","Apr","May","Jun")
starlabs = rep(starlab, 8)
yr = c("`05", "`06", "`07", "`08", "`09", "`10", "`11")
yrs = c("`04","`04","`04","`04","`04","`04", rep(yr, each = 12), "`12","`12","`12","`12","`12","`12")
starlabsyr = paste(starlabs, yrs, sep = '')

#export table to latex format
rownames(stargaze) = starlabsyr
colnames(stargaze) = c('$EMA_k$', 'SE')
print(xtable(stargaze), include.rownames = T, sanitize.colnames.function =identity)

#----------------------------------------------EMA Plot-----------------------------------------------------------# 
#
# Create the EMA plot (Figure 2 in paper)
#
#------------------------------------------------------------------------------------------------------------------#

month = c("Jul","Aug","Sep","Oct","Nov","Dec", "Jan","Feb","Mar","Apr","May","Jun")
month <- rep(month, 8)
ema.ts <- data.frame(ema.beta, month)
date <- as.yearmon(2004 + seq(6,101)/12)
ema.ts$date <- as.Date(date)

emaplot1 <- ggplot(ema.ts,aes(date, ema.beta, label=date)) +
  geom_point(colour="black", size = 4, pch = 20) +
  geom_point(data=subset(ema.ts, month == "Dec"), aes(date, ema.beta, label = date), colour = "red", size = 4, pch =20) +
  geom_point(data=ema.ts[54,], aes(date, ema.beta, label = date), colour="red", size = 5, pch = 17) +
  geom_point(data=subset(ema.ts, month == "Apr"), aes(date, ema.beta, label = date), colour = "blue", size = 4, pch =20) +
  geom_point(data=ema.ts[58,], aes(date, ema.beta, label = date), colour="blue", size = 5, pch = 17) +
  geom_point(data=subset(ema.ts, month == "May"), aes(date, ema.beta, label = date), colour = "orange", size = 4, pch =20) +
  geom_point(data=ema.ts[59,], aes(date, ema.beta, label = date), colour="orange", size = 5, pch = 17) +
  geom_point(data=ema.ts[95,], aes(date, ema.beta, label = date), colour="orange", size = 5, pch = 15) +
  geom_point(data=ema.ts[96,], aes(date, ema.beta, label = date), colour="black", size = 5, pch = 15) +
  geom_point(data=subset(ema.ts, month == "Mar"), aes(date, ema.beta, label = date), colour = "purple", size = 4, pch =20) +
  geom_point(data=ema.ts[57,], aes(date, ema.beta, label = date), colour="purple", size = 5, pch = 17) +
  xlab("") +
  ylab("") +
  theme(axis.title.y = element_text(family="serif", face="bold", size=14)) +
  geom_text(data=ema.ts[54, ], label="December '08", family="serif", face="bold", size=5, vjust=.35, hjust = -.15, colour = "black") +
  geom_text(data=ema.ts[57, ], label="March '09", family="serif", face="bold", size=5, vjust=.35, hjust = -.15, colour = "black") +
  geom_text(data=ema.ts[58, ], label="April '09", family="serif", face="bold", size=5, vjust=.35, hjust = -.25, colour = "black") +
  geom_text(data=ema.ts[59, ], label="May '09", family="serif", face="bold", size=5, vjust=.35, hjust = -.25, colour = "black") +
  theme(panel.background = element_blank()) +
  scale_x_date(breaks = "1 year", labels=date_format("%Y")) +
  scale_y_continuous(breaks=seq(-500,900, 200)) +

  theme(axis.line = element_line(colour = "black"),  
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank()) 

print(emaplot1)

ggsave(plot=emaplot1, filename=paste('emaplotcolorV3', '.pdf', sep=''), scale=1.2,height=5, width=7)



#------------------------------EMA*LGA regression----------------------------------------------#
#
# This is the monster regression used in the paper. It generates 5952 (96*62) coefficients. 
# These coefficients are used to make the density plots in the paper (Figures 3 and 6)
#
#----------------------------------------------------------------------------------------------#

#create 96 EMA variables to match the dummies in the dataset 
e = rep(1:96)
emalga_variables = paste('EMA_', e, sep = "")

#original regression formula: real monthly expenditures per egm on monthly dummies using 'Jan' as 
# the base + a linear time trend + LGA fixed effects 
emalga_regression <- rMonExp_EGM~ Month2+Month3+Month4+Month5+Month6+Month7+Month8+Month9+
  Month10+Month11+Month12+Yrs_minus_2004 

#Create list of 96 copies of the original regression each w/ a different EMA dummy 
emalga_regression_list <- lapply(emalga_variables, function(x, orig = emalga_regression) { 
  new <- reformulate(c( x, '.'))
  update(orig, new)})

#apply 96 models to each LGA by subsetting the data, each model has a dummy for t=1:96
Alpine               <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short1  == 1,])
Ararat               <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short2  == 1,])
Ballarat             <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short3  == 1,])
Banyule              <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short4  == 1,])
Bass_Coast           <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short5  == 1,])
Baw_Baw              <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short6  == 1,])
Bayside              <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short7  == 1,])
Benalla              <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short8  == 1,])
Boroondara           <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short9  == 1,])
Brimbank             <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short10 == 1,])
Campaspe             <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short11 == 1,])
Cardinia             <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short12 == 1,])
Casey                <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short13 == 1,])
Central_Goldfields   <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short14 == 1,])
Colac_Otway          <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short15 == 1,])
Corangamite          <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short16 == 1,])
Darebin              <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short17 == 1,])
East_Gippsland       <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short18 == 1,])
Frankston            <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short19 == 1,])
Glen_Eira            <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short20 == 1,])
Glenelg              <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short21 == 1,])
Greater_Bendigo      <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short22 == 1,])
Greater_Dandenong    <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short23 == 1,])
Greater_Geelong      <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short24 == 1,])
Greater_Shepparton   <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short25 == 1,])
Hobsons_Bay          <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short26 == 1,])
Horsham              <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short27 == 1,])
Hume                 <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short28 == 1,])
Kingston             <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short29 == 1,])
Knox                 <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short30 == 1,])
Latrobe              <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short31 == 1,])
Macedon_Ranges       <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short32 == 1,])
Manningham           <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short33 == 1,])
Mansfield            <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short34 == 1,])
Maribyrnong          <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short35 == 1,])
Maroondah            <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short36 == 1,])
Melbourne            <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short37 == 1,])
Melton               <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short38 == 1,])
Mildura              <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short39 == 1,])
Mitchell             <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short40 == 1,])
Monash               <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short41 == 1,])
Moonee_Valley        <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short42 == 1,])
Moorabool            <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short43 == 1,])
Moreland             <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short44 == 1,])
Mornington_Peninsula <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short45 == 1,])
Murrindindi          <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short46 == 1,])
Nillumbik            <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short47 == 1,])
Northern_Grampians   <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short48 == 1,])
Port_Phillip         <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short49 == 1,])
South_Gippsland      <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short50 == 1,])
Stonnington          <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short51 == 1,])
Surf_Coast           <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short52 == 1,])
Swan_Hill            <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short53 == 1,])
Wangaratta           <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short54 == 1,])
Warrnambool          <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short55 == 1,])
Wellington           <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short56 == 1,])
Whitehorse           <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short57 == 1,])
Whittlesea           <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short58 == 1,])
Wodonga              <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short59 == 1,])
Wyndham              <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short60 == 1,])
Yarra                <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short61 == 1,])
Yarra_Ranges         <- lapply(emalga_regression_list, lm, data = pokies[pokies$LGA_short62 == 1,])

#name the models to correspond to each LGA; each LGA has 96 regressions attached to it as per 
# the object 'emalga_regression_list' above. 
models = c(Alpine,Ararat,Ballarat,Banyule,Bass_Coast,Baw_Baw,Bayside,Benalla,Boroondara,           
    Brimbank,Campaspe,Cardinia,Casey,Central_Goldfields,Colac_Otway,Corangamite,         
    Darebin,East_Gippsland,Frankston,Glen_Eira,Glenelg,Greater_Bendigo,
    Greater_Dandenong,Greater_Geelong,Greater_Shepparton,Hobsons_Bay,Horsham,Hume,                 
    Kingston,Knox,Latrobe,Macedon_Ranges,Manningham,Mansfield,Maribyrnong,Maroondah,            
    Melbourne,Melton,Mildura,Mitchell,Monash,Moonee_Valley,Moorabool,Moreland,             
    Mornington_Peninsula,Murrindindi,Nillumbik,Northern_Grampians,Port_Phillip,         
    South_Gippsland,Stonnington,Surf_Coast,Swan_Hill,Wangaratta,Warrnambool,          
    Wellington,Whitehorse,Whittlesea,Wodonga,Wyndham,Yarra,Yarra_Ranges)

hnam = c('Alpine','Ararat','Ballarat','Banyule','Bass_Coast','Baw_Baw','Bayside','Benalla','Boroondara',           
'Brimbank','Campaspe','Cardinia','Casey','Central_Goldfields','Colac_Otway','Corangamite',         
'Darebin','East_Gippsland','Frankston','Glen_Eira','Glenelg','Greater_Bendigo',
'Greater_Dandenong','Greater_Geelong','Greater_Shepparton','Hobsons_Bay','Horsham','Hume',                 
'Kingston','Knox','Latrobe','Macedon_Ranges','Manningham','Mansfield','Maribyrnong','Maroondah',            
'Melbourne','Melton','Mildura','Mitchell','Monash','Moonee_Valley','Moorabool','Moreland',             
'Mornington_Peninsula','Murrindindi','Nillumbik','Northern_Grampians','Port_Phillip',         
'South_Gippsland','Stonnington','Surf_Coast','Swan_Hill','Wangaratta','Warrnambool',          
'Wellington','Whitehorse','Whittlesea','Wodonga','Wyndham','Yarra','Yarra_Ranges')

mod = rep(hnam, each =96)
modnames = paste(mod, '_period', 1:96, sep="")
names(models) = modnames 

#create summarries for all the models
summaries <- lapply(models, summary) 

#Extract the coefficients from the models 
coefficients <- lapply(models, coef)

#Create table with coefficient estimates and standard errors for all 96 models
coef_tables <- lapply(summaries, '[[', 'coefficients')

#Extract 5952 coefficients 
emalga.beta <- matrix(nrow=5952,ncol=1) 
names(models) <-paste('reg.out', 1:5952, sep="")

for (i in 1:5952){
  emalga.beta[i,] <- c(models[[paste('reg.out',i,sep='')]]$coefficients[2])  
}

#a nice property, the mean is approx. zero! phew! 
mean(emalga.beta)

#rows are alpine, 1:96 , Ararat 1:96 etc. ; generate sequence to identify
#monthly EMAs by LGA (i.e. Dec EMAs start at 6 and increase by 12 until)
#March, April, May, June, December and August are the odd ones 
juls = seq(1, 5952, by = 12)
augs = seq(2, 5952, by = 12)
seps = seq(3, 5952, by = 12)
octs = seq(4, 5952, by = 12)
novs = seq(5, 5952, by = 12)
decs = seq(6, 5952, by = 12)
jans = seq(7, 5952, by = 12)
febs = seq(8, 5952, by = 12)
mars = seq(9, 5952, by = 12)
aprs = seq(10, 5952, by = 12)
mays = seq(11, 5952, by = 12)
juns = seq(12, 5952, by = 12)
  

#Extract EMAs 
jul_emas = as.matrix(emalga.beta[juls,1]) 
aug_emas = as.matrix(emalga.beta[augs,1]) 
sep_emas = as.matrix(emalga.beta[seps,1]) 
oct_emas = as.matrix(emalga.beta[octs,1]) 
nov_emas = as.matrix(emalga.beta[novs,1]) 
dec_emas = as.matrix(emalga.beta[decs,1]) 
jan_emas = as.matrix(emalga.beta[jans,1])
feb_emas = as.matrix(emalga.beta[febs,1])
mar_emas = as.matrix(emalga.beta[mars,1])
apr_emas = as.matrix(emalga.beta[aprs,1])
may_emas = as.matrix(emalga.beta[mays,1])
jun_emas = as.matrix(emalga.beta[juns,1])


#----------------------------------Plots for each month--------------------------------------------------#
#
# This section creates density plots of the coefficients for each month in the sample. That is, we have 
# 96 density curves, each estimated using 62 datapoints (one coefficient from each LGA). These plots are
# used to produce Figures 3 and 6 in the paper. 
#
#--------------------------------------------------------------------------------------------------------#



#--------------------------------------------December----------------------------------------------------#
#create indicators for Dec04 ... Dec11, then attach to data frame w/ estimates
declab = c('Dec04', 'Dec05', 'Dec06', 'Dec07', 'Dec08', 'Dec09', 'Dec10', 'Dec11')
declabs = rep(declab, 62)

rownames(dec_emas) = declabs
colnames(dec_emas) = 'EMA'

#add in factor ID for the 8 levels 
dec04 = as.numeric(rownames(dec_emas) == 'Dec04')
dec05 = as.numeric(rownames(dec_emas) == 'Dec05')
dec06 = as.numeric(rownames(dec_emas) == 'Dec06')
dec07 = as.numeric(rownames(dec_emas) == 'Dec07')
dec08 = as.numeric(rownames(dec_emas) == 'Dec08')
dec09 = as.numeric(rownames(dec_emas) == 'Dec09')
dec10 = as.numeric(rownames(dec_emas) == 'Dec10')
dec11 = as.numeric(rownames(dec_emas) == 'Dec11')

dec_ema = as.numeric(dec_emas)
dec_ema_df = data.frame(dec_ema, dec04,dec05,dec06,dec07,dec08,dec09,dec10,dec11)

dec_ema_df$new[dec_ema_df$dec04==1]<-4
dec_ema_df$new[dec_ema_df$dec05==1]<-5
dec_ema_df$new[dec_ema_df$dec06==1]<-6
dec_ema_df$new[dec_ema_df$dec07==1]<-7
dec_ema_df$new[dec_ema_df$dec08==1]<-8
dec_ema_df$new[dec_ema_df$dec09==1]<-9
dec_ema_df$new[dec_ema_df$dec10==1]<-10
dec_ema_df$new[dec_ema_df$dec11==1]<-11


#using GGPLOT2... WORD! 
library(ggplot2)
library('RColorBrewer')


#December density plot  
cbbPalette <- brewer.pal(8, "Accent")


dec.dens1 = ggplot(dec_ema_df, aes(x=dec_ema, group = new)) + 
  geom_density(aes(colour=factor(new), linetype = factor(new)), size = 1, alpha = .60, kernel='epanechnikov') + 
  xlab(NULL) + ylab(NULL) + 
  scale_linetype_manual(values = c('solid', 'solid','solid','solid','longdash','solid','solid','solid'),  
                        labels=c('`04','`05','`06','`07','`08','`09','`10','`11')) +
                          scale_colour_manual(values = cbbPalette,labels= c('`04','`05','`06','`07','`08','`09','`10','`11'),name= NULL) +            
                          scale_y_continuous(limits=c(0.0000, 0.0030), labels = percent) +
                          scale_x_continuous(limits=c(-2000,2000)) +
                          theme(axis.line = element_line(colour = "black"),
                                panel.grid.major = element_blank(),
                                panel.grid.minor = element_blank(),
                                panel.border = element_blank(),
                                panel.background = element_blank(),
                                legend.title=element_blank(),
                                legend.position = 'bottom', 
                                legend.text = element_text(color= 'black', size = 12, face = 'bold'),
                                legend.key = element_rect(fill = 'white', size = 0.5, linetype='solid'),
                                legend.key.size = unit(1.5, 'lines'))


print(dec.dens1)

setwd('c:/Users/peytonk/Dropbox/Pokies/Replication')
ggsave(plot=dec.dens1, filename=paste('decdensV2', '.pdf', sep=''), scale=1,height=5, width=7)


#--------------------------------------------January----------------------------------------------------------#

janlab = c('Jan05', 'Jan06', 'Jan07', 'Jan08', 'Jan09', 'Jan10', 'Jan11', 'Jan12')
janlabs = rep(janlab, 62)

rownames(jan_emas) = janlabs   
colnames(jan_emas) = 'EMA'

#add in factor ID for the 8 levels 
jan05 = as.numeric(rownames(jan_emas) == 'Jan05')
jan06 = as.numeric(rownames(jan_emas) == 'Jan06')
jan07 = as.numeric(rownames(jan_emas) == 'Jan07')
jan08 = as.numeric(rownames(jan_emas) == 'Jan08')
jan09 = as.numeric(rownames(jan_emas) == 'Jan09')
jan10 = as.numeric(rownames(jan_emas) == 'Jan10')
jan11 = as.numeric(rownames(jan_emas) == 'Jan11')
jan12 = as.numeric(rownames(jan_emas) == 'Jan12')

jan_ema = as.numeric(jan_emas)
jan_ema_df = data.frame(jan_ema, jan05,jan06,jan07,jan08,jan09,jan10,jan11,jan12)

jan_ema_df$new[jan_ema_df$jan05==1]<-5
jan_ema_df$new[jan_ema_df$jan06==1]<-6
jan_ema_df$new[jan_ema_df$jan07==1]<-7
jan_ema_df$new[jan_ema_df$jan08==1]<-8
jan_ema_df$new[jan_ema_df$jan09==1]<-9
jan_ema_df$new[jan_ema_df$jan10==1]<-10
jan_ema_df$new[jan_ema_df$jan11==1]<-11
jan_ema_df$new[jan_ema_df$jan12==1]<-12


jan.dens1 = ggplot(jan_ema_df, aes(x=jan_ema, group = new)) + 
  geom_density(aes(colour=factor(new)), size = 1, alpha = .70,
               kernel='epanechnikov') + 
                 xlab(NULL) + ylab(NULL) + 
                 scale_colour_manual(values = cbbPalette,
                                     labels= c('`05','`06','`07','`08','`09','`10','`11', '`12'),
                                     name= NULL) +
                                       scale_y_continuous(limits=c(0.0000, 0.0025), labels = percent) +
                                       scale_x_continuous(limits=c(-2000,2000)) +
                                       theme(axis.line = element_line(colour = "black"),
                                             panel.grid.major = element_blank(),
                                             panel.grid.minor = element_blank(),
                                             panel.border = element_blank(),
                                             panel.background = element_blank(),
                                             legend.title=element_blank(),
                                             legend.position = 'bottom', 
                                             legend.text = element_text(color= 'black', size = 12, face = 'bold'),
                                             legend.key = element_rect(fill = 'white', size = 0.5, linetype='solid'),
                                             legend.key.size = unit(1.5, 'lines'))

print(jan.dens1)


setwd('c:/Users/peytonk/Dropbox/Pokies/Replication')
ggsave(plot=jan.dens1, filename=paste('jandensV2', '.pdf', sep=''), scale=1,height=5, width=7)


#--------------------------------------------------February--------------------------------------------------------#

feblab = c('Feb05', 'Feb06', 'Feb07', 'Feb08', 'Feb09', 'Feb10', 'Feb11', 'Feb12')
feblabs = rep(feblab, 62)

rownames(feb_emas) = feblabs   
colnames(feb_emas) = 'EMA'

#add in factor ID for the 8 levels 
feb05 = as.numeric(rownames(feb_emas) == 'Feb05')
feb06 = as.numeric(rownames(feb_emas) == 'Feb06')
feb07 = as.numeric(rownames(feb_emas) == 'Feb07')
feb08 = as.numeric(rownames(feb_emas) == 'Feb08')
feb09 = as.numeric(rownames(feb_emas) == 'Feb09')
feb10 = as.numeric(rownames(feb_emas) == 'Feb10')
feb11 = as.numeric(rownames(feb_emas) == 'Feb11')
feb12 = as.numeric(rownames(feb_emas) == 'Feb12')

feb_ema = as.numeric(feb_emas)
feb_ema_df = data.frame(feb_ema, feb05,feb06,feb07,feb08,feb09,feb10,feb11,feb12)

feb_ema_df$new[feb_ema_df$feb05==1]<-5
feb_ema_df$new[feb_ema_df$feb06==1]<-6
feb_ema_df$new[feb_ema_df$feb07==1]<-7
feb_ema_df$new[feb_ema_df$feb08==1]<-8
feb_ema_df$new[feb_ema_df$feb09==1]<-9
feb_ema_df$new[feb_ema_df$feb10==1]<-10
feb_ema_df$new[feb_ema_df$feb11==1]<-11
feb_ema_df$new[feb_ema_df$feb12==1]<-12



feb.dens1 = ggplot(feb_ema_df, aes(x=feb_ema, group = new)) + 
  geom_density(aes(colour=factor(feb_ema_df$new)), size = 1, alpha = .70,kernel='epanechnikov') +
  xlab(NULL) + ylab(NULL) + 
  scale_colour_manual(values = cbbPalette,
                      labels= c('`05','`06','`07','`08','`09','`10','`11', '`12'),
                      name= NULL) +
                        scale_y_continuous(limits=c(0.0000, 0.0025), labels = percent) +
                        scale_x_continuous(limits=c(-2000,2000)) +
                        theme(axis.line = element_line(colour = "black"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.border = element_blank(),
                              panel.background = element_blank(),
                              legend.title=element_blank(),
                              legend.position = 'bottom', 
                              legend.text = element_text(color= 'black', size = 12, face = 'bold'),
                              legend.key = element_rect(fill = 'white', size = 0.5, linetype='solid'),
                              legend.key.size = unit(1.5, 'lines'))

print(feb.dens1)


setwd('c:/Users/peytonk/Dropbox/Pokies/Replication')
ggsave(plot=feb.dens1, filename=paste('febdensV2', '.pdf', sep=''), scale=1,height=5, width=7)




#-------------------------------------------------March---------------------------------------------------------#

marlab = c('Mar05', 'Mar06', 'Mar07', 'Mar08', 'Mar09', 'Mar10', 'Mar11', 'Mar12')
marlabs = rep(marlab, 62)

rownames(mar_emas) = marlabs   
colnames(mar_emas) = 'EMA'

#add in factor ID for the 8 levels 
mar05 = as.numeric(rownames(mar_emas) == 'Mar05')
mar06 = as.numeric(rownames(mar_emas) == 'Mar06')
mar07 = as.numeric(rownames(mar_emas) == 'Mar07')
mar08 = as.numeric(rownames(mar_emas) == 'Mar08')
mar09 = as.numeric(rownames(mar_emas) == 'Mar09')
mar10 = as.numeric(rownames(mar_emas) == 'Mar10')
mar11 = as.numeric(rownames(mar_emas) == 'Mar11')
mar12 = as.numeric(rownames(mar_emas) == 'Mar12')

mar_ema = as.numeric(mar_emas)
mar_ema_df = data.frame(mar_ema, mar05,mar06,mar07,mar08,mar09,mar10,mar11,mar12)

mar_ema_df$new[mar_ema_df$mar05==1]<-5
mar_ema_df$new[mar_ema_df$mar06==1]<-6
mar_ema_df$new[mar_ema_df$mar07==1]<-7
mar_ema_df$new[mar_ema_df$mar08==1]<-8
mar_ema_df$new[mar_ema_df$mar09==1]<-9
mar_ema_df$new[mar_ema_df$mar10==1]<-10
mar_ema_df$new[mar_ema_df$mar11==1]<-11
mar_ema_df$new[mar_ema_df$mar12==1]<-12


mar.dens1 = ggplot(mar_ema_df, aes(x=mar_ema, group = new)) + 
  geom_density(aes(colour=factor(new), linetype= factor(new)), size = 1, alpha = .70, kernel = 'epanechnikov') + 
  xlab(NULL) + ylab(NULL) + 
  scale_linetype_manual(values=c('solid', 'solid','solid','solid','longdash','solid','solid','solid'),
                        labels=c('`05','`06','`07','`08','`09','`10','`11','`12')) +
                          scale_colour_manual(values = cbbPalette,labels=c('`05','`06','`07','`08','`09','`10','`11','`12'), name=NULL) +
                          scale_y_continuous(limits=c(0.0000, 0.0030), labels = percent) +
                          scale_x_continuous(limits=c(-2000,2000)) +
                          theme(axis.line = element_line(colour = "black"),
                                panel.grid.major = element_blank(),
                                panel.grid.minor = element_blank(),
                                panel.border = element_blank(),
                                panel.background = element_blank(),
                                legend.title=element_blank(),
                                legend.position = 'bottom', 
                                legend.text = element_text(color= 'black', size = 12, face = 'bold'),
                                legend.key = element_rect(fill = 'white', size = 0.5, linetype='solid'),
                                legend.key.size = unit(1.5, 'lines'))

print(mar.dens1)


setwd('c:/Users/peytonk/Dropbox/Pokies/Replication')
ggsave(plot=mar.dens1, filename=paste('mardensV2', '.pdf', sep=''), scale=1,height=5, width=7)


#-------------------------------------------------April----------------------------------------------------------#


aprlab = c('Apr05', 'Apr06', 'Apr07', 'Apr08', 'Apr09', 'Apr10', 'Apr11', 'Apr12')
aprlabs = rep(aprlab, 62)

rownames(apr_emas) = aprlabs   
colnames(apr_emas) = 'EMA'

#add in factor ID for the 8 levels 
apr05 = as.numeric(rownames(apr_emas) == 'Apr05')
apr06 = as.numeric(rownames(apr_emas) == 'Apr06')
apr07 = as.numeric(rownames(apr_emas) == 'Apr07')
apr08 = as.numeric(rownames(apr_emas) == 'Apr08')
apr09 = as.numeric(rownames(apr_emas) == 'Apr09')
apr10 = as.numeric(rownames(apr_emas) == 'Apr10')
apr11 = as.numeric(rownames(apr_emas) == 'Apr11')
apr12 = as.numeric(rownames(apr_emas) == 'Apr12')

apr_ema = as.numeric(apr_emas)
apr_ema_df = data.frame(apr_ema, apr05,apr06,apr07,apr08,apr09,apr10,apr11,apr12)

apr_ema_df$new[apr_ema_df$apr05==1]<-5
apr_ema_df$new[apr_ema_df$apr06==1]<-6
apr_ema_df$new[apr_ema_df$apr07==1]<-7
apr_ema_df$new[apr_ema_df$apr08==1]<-8
apr_ema_df$new[apr_ema_df$apr09==1]<-9
apr_ema_df$new[apr_ema_df$apr10==1]<-10
apr_ema_df$new[apr_ema_df$apr11==1]<-11
apr_ema_df$new[apr_ema_df$apr12==1]<-12



apr.dens1 = ggplot(apr_ema_df, aes(x=apr_ema, group = new)) + 
  geom_density(aes(colour=factor(new), linetype= factor(new)), size = 1, alpha = .70, kernel='epanechnikov') + 
  xlab(NULL) + ylab(NULL) + 
  scale_linetype_manual(values=c('solid', 'solid','solid','solid','longdash','solid','solid','solid'),
                        labels=c('`05','`06','`07','`08','`09','`10','`11','`12')) +
                          scale_colour_manual(values = cbbPalette,labels=c('`05','`06','`07','`08','`09','`10','`11','`12'), name=NULL) +
                          scale_y_continuous(limits=c(0.0000, 0.0030), labels = percent) +
                          scale_x_continuous(limits=c(-2000,2000)) +
                          theme(axis.line = element_line(colour = "black"),
                                panel.grid.major = element_blank(),
                                panel.grid.minor = element_blank(),
                                panel.border = element_blank(),
                                panel.background = element_blank(),
                                legend.title=element_blank(),
                                legend.position = 'bottom', 
                                legend.text = element_text(color= 'black', size = 12, face = 'bold'),
                                legend.key = element_rect(fill = 'white', size = 0.5, linetype='solid'),
                                legend.key.size = unit(1.5, 'lines'))

print(apr.dens1)



setwd('c:/Users/peytonk/Dropbox/Pokies/Replication')
ggsave(plot=apr.dens1, filename=paste('aprdensV2', '.pdf', sep=''), scale=1,height=5, width=7)


#----------------------------------------------------May-------------------------------------------------------#


maylab = c('May05', 'May06', 'May07', 'May08', 'May09', 'May10', 'May11', 'May12')
maylabs = rep(maylab, 62)

rownames(may_emas) = maylabs   
colnames(may_emas) = 'EMA'

#add in factor ID for the 8 levels 
may05 = as.numeric(rownames(may_emas) == 'May05')
may06 = as.numeric(rownames(may_emas) == 'May06')
may07 = as.numeric(rownames(may_emas) == 'May07')
may08 = as.numeric(rownames(may_emas) == 'May08')
may09 = as.numeric(rownames(may_emas) == 'May09')
may10 = as.numeric(rownames(may_emas) == 'May10')
may11 = as.numeric(rownames(may_emas) == 'May11')
may12 = as.numeric(rownames(may_emas) == 'May12')

may_ema = as.numeric(may_emas)
may_ema_df = data.frame(may_ema, may05,may06,may07,may08,may09,may10,may11,may12)

may_ema_df$new[may_ema_df$may05==1]<-5
may_ema_df$new[may_ema_df$may06==1]<-6
may_ema_df$new[may_ema_df$may07==1]<-7
may_ema_df$new[may_ema_df$may08==1]<-8
may_ema_df$new[may_ema_df$may09==1]<-9
may_ema_df$new[may_ema_df$may10==1]<-10
may_ema_df$new[may_ema_df$may11==1]<-11
may_ema_df$new[may_ema_df$may12==1]<-12



may.dens1 = ggplot(may_ema_df, aes(x=may_ema, group = new)) + 
  geom_density(aes(colour=factor(new), linetype = factor(new)), size = 1, alpha = .70,kernel='epanechnikov') + 
  xlab(NULL) + ylab(NULL) + 
  scale_linetype_manual(values=c('solid', 'solid','solid','solid','longdash','solid','solid','longdash'),
                        labels=c('`05','`06','`07','`08','`09','`10','`11','`12')) +
                          scale_colour_manual(values = cbbPalette,labels=c('`05','`06','`07','`08','`09','`10','`11','`12'), name=NULL) +
                          scale_y_continuous(limits=c(0.0000, 0.0030), labels = percent) +
                          scale_x_continuous(limits=c(-2000,2000)) +
                          theme(axis.line = element_line(colour = "black"),
                                panel.grid.major = element_blank(),
                                panel.grid.minor = element_blank(),
                                panel.border = element_blank(),
                                panel.background = element_blank(),
                                legend.title=element_blank(),
                                legend.position = 'bottom', 
                                legend.text = element_text(color= 'black', size = 12, face = 'bold'),
                                legend.key = element_rect(fill = 'white', size = 0.5, linetype='solid'),
                                legend.key.size = unit(1.5, 'lines'))

print(may.dens1)



setwd('c:/Users/peytonk/Dropbox/Pokies/Replication')
ggsave(plot=may.dens1, filename=paste('maydensV2', '.pdf', sep=''), scale=1,height=5, width=7)


#------------------------------------------------------June---------------------------------------------------------#


junlab = c('Jun05', 'Jun06', 'Jun07', 'Jun08', 'Jun09', 'Jun10', 'Jun11', 'Jun12')
junlabs = rep(junlab, 62)

rownames(jun_emas) = junlabs   
colnames(jun_emas) = 'EMA'

#add in factor ID for the 8 levels 
jun05 = as.numeric(rownames(jun_emas) == 'Jun05')
jun06 = as.numeric(rownames(jun_emas) == 'Jun06')
jun07 = as.numeric(rownames(jun_emas) == 'Jun07')
jun08 = as.numeric(rownames(jun_emas) == 'Jun08')
jun09 = as.numeric(rownames(jun_emas) == 'Jun09')
jun10 = as.numeric(rownames(jun_emas) == 'Jun10')
jun11 = as.numeric(rownames(jun_emas) == 'Jun11')
jun12 = as.numeric(rownames(jun_emas) == 'Jun12')

jun_ema = as.numeric(jun_emas)
jun_ema_df = data.frame(jun_ema, jun05,jun06,jun07,jun08,jun09,jun10,jun11,jun12)

jun_ema_df$new[jun_ema_df$jun05==1]<-5
jun_ema_df$new[jun_ema_df$jun06==1]<-6
jun_ema_df$new[jun_ema_df$jun07==1]<-7
jun_ema_df$new[jun_ema_df$jun08==1]<-8
jun_ema_df$new[jun_ema_df$jun09==1]<-9
jun_ema_df$new[jun_ema_df$jun10==1]<-10
jun_ema_df$new[jun_ema_df$jun11==1]<-11
jun_ema_df$new[jun_ema_df$jun12==1]<-12



jun.dens1 = ggplot(jun_ema_df, aes(x=jun_ema, group = new)) + 
  geom_density(aes(colour=factor(new), linetype = factor(new)), size = 1, alpha = .70,kernel='epanechnikov') + 
  xlab(NULL) + ylab(NULL) + 
  scale_linetype_manual(values=c('solid', 'solid','solid','solid','longdash','solid','solid','longdash'),
                        labels=c('`05','`06','`07','`08','`09','`10','`11','`12')) +
                          scale_colour_manual(values = cbbPalette,labels=c('`05','`06','`07','`08','`09','`10','`11','`12'), name=NULL) +
                          scale_y_continuous(limits=c(0.0000, 0.0030), labels = percent) +
                          scale_x_continuous(limits=c(-2000,2000)) +
                          theme(axis.line = element_line(colour = "black"),
                                panel.grid.major = element_blank(),
                                panel.grid.minor = element_blank(),
                                panel.border = element_blank(),
                                panel.background = element_blank(),
                                legend.title=element_blank(),
                                legend.position = 'bottom', 
                                legend.text = element_text(color= 'black', size = 12, face = 'bold'),
                                legend.key = element_rect(fill = 'white', size = 0.5, linetype='solid'),
                                legend.key.size = unit(1.5, 'lines')) 


print(jun.dens1)


setwd('c:/Users/peytonk/Dropbox/Pokies/Replication')
ggsave(plot=jun.dens1, filename=paste('jundensV2', '.pdf', sep=''), scale=1,height=5, width=7)


#-----------------------------------------------------July---------------------------------------------------------#

#create indicators for Jul04 ... Jul11, then attach to data frame w/ estimates
jullab = c('Jul04', 'Jul05', 'Jul06', 'Jul07', 'Jul08', 'Jul09', 'Jul10', 'Jul11')
jullabs = rep(jullab, 62)

rownames(jul_emas) = jullabs
colnames(jul_emas) = 'EMA'

#add in factor ID for the 8 levels 
jul04 = as.numeric(rownames(jul_emas) == 'Jul04')
jul05 = as.numeric(rownames(jul_emas) == 'Jul05')
jul06 = as.numeric(rownames(jul_emas) == 'Jul06')
jul07 = as.numeric(rownames(jul_emas) == 'Jul07')
jul08 = as.numeric(rownames(jul_emas) == 'Jul08')
jul09 = as.numeric(rownames(jul_emas) == 'Jul09')
jul10 = as.numeric(rownames(jul_emas) == 'Jul10')
jul11 = as.numeric(rownames(jul_emas) == 'Jul11')

jul_ema = as.numeric(jul_emas)
jul_ema_df = data.frame(jul_ema, jul04,jul05,jul06,jul07,jul08,jul09,jul10,jul11)

jul_ema_df$new[jul_ema_df$jul04==1]<-4
jul_ema_df$new[jul_ema_df$jul05==1]<-5
jul_ema_df$new[jul_ema_df$jul06==1]<-6
jul_ema_df$new[jul_ema_df$jul07==1]<-7
jul_ema_df$new[jul_ema_df$jul08==1]<-8
jul_ema_df$new[jul_ema_df$jul09==1]<-9
jul_ema_df$new[jul_ema_df$jul10==1]<-10
jul_ema_df$new[jul_ema_df$jul11==1]<-11

jul.dens1 = ggplot(jul_ema_df, aes(x=jul_ema, group = new)) + 
  geom_density(aes(colour=factor(new)), size = 1, alpha = .60,kernel='epanechnikov') +
  xlab(NULL) + ylab(NULL) + 
  scale_colour_manual(values = cbbPalette,
                      labels= c('`04','`05','`06','`07','`08','`09','`10','`11'),
                      name= '') +
                        scale_y_continuous(limits=c(0.0000, 0.0030), labels = percent) +
                        scale_x_continuous(limits=c(-2000,2000)) +
                        theme(axis.line = element_line(colour = "black"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.border = element_blank(),
                              panel.background = element_blank(),
                              legend.title=element_blank(),
                              legend.position = 'bottom', 
                              legend.text = element_text(color= 'black', size = 12, face = 'bold'),
                              legend.key = element_rect(fill = 'white', size = 0.5, linetype='solid'),
                              legend.key.size = unit(1.5, 'lines'))


print(jul.dens1)


setwd('c:/Users/peytonk/Dropbox/Pokies/Replication')
ggsave(plot=jul.dens1, filename=paste('juldensV2', '.pdf', sep=''), scale=1,height=5, width=7)


#----------------------------------------------------August--------------------------------------------------------#

#create indicators for Aug04 ... Aug11, then attach to data frame w/ estimates
auglab = c('Aug04', 'Aug05', 'Aug06', 'Aug07', 'Aug08', 'Aug09', 'Aug10', 'Aug11')
auglabs = rep(auglab, 62)

rownames(aug_emas) = auglabs
colnames(aug_emas) = 'EMA'

#add in factor ID for the 8 levels 
aug04 = as.numeric(rownames(aug_emas) == 'Aug04')
aug05 = as.numeric(rownames(aug_emas) == 'Aug05')
aug06 = as.numeric(rownames(aug_emas) == 'Aug06')
aug07 = as.numeric(rownames(aug_emas) == 'Aug07')
aug08 = as.numeric(rownames(aug_emas) == 'Aug08')
aug09 = as.numeric(rownames(aug_emas) == 'Aug09')
aug10 = as.numeric(rownames(aug_emas) == 'Aug10')
aug11 = as.numeric(rownames(aug_emas) == 'Aug11')

aug_ema = as.numeric(aug_emas)
aug_ema_df = data.frame(aug_ema, aug04,aug05,aug06,aug07,aug08,aug09,aug10,aug11)

aug_ema_df$new[aug_ema_df$aug04==1]<-4
aug_ema_df$new[aug_ema_df$aug05==1]<-5
aug_ema_df$new[aug_ema_df$aug06==1]<-6
aug_ema_df$new[aug_ema_df$aug07==1]<-7
aug_ema_df$new[aug_ema_df$aug08==1]<-8
aug_ema_df$new[aug_ema_df$aug09==1]<-9
aug_ema_df$new[aug_ema_df$aug10==1]<-10
aug_ema_df$new[aug_ema_df$aug11==1]<-11


#August density plot  
aug.dens1 = ggplot(aug_ema_df, aes(x=aug_ema, group = new)) + 
  geom_density(aes(colour=factor(new)), size = 1, alpha = .60,kernel='epanechnikov') +
  xlab(NULL) + ylab(NULL) +
  scale_colour_manual(values = cbbPalette,
                      labels= c('`04','`05','`06','`07','`08','`09','`10','`11'),
                      name= '') +
                        scale_y_continuous(limits=c(0.0000, 0.0030), labels = percent) +
                        scale_x_continuous(limits=c(-2000,2000)) +
                        theme(axis.line = element_line(colour = "black"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.border = element_blank(),
                              panel.background = element_blank(),
                              legend.title=element_blank(),
                              legend.position = 'bottom', 
                              legend.text = element_text(color= 'black', size = 12, face = 'bold'),
                              legend.key = element_rect(fill = 'white', size = 0.5, linetype='solid'),
                              legend.key.size = unit(1.5, 'lines'))


print(aug.dens1)

setwd('c:/Users/peytonk/Dropbox/Pokies/Replication')
ggsave(plot=aug.dens1, filename=paste('augdensV2', '.pdf', sep=''), scale=1,height=5, width=7)



#-------------------------------------------------September-------------------------------------------------------#
#create indicators for Sep04 ... Sep11, then attach to data frame w/ estimates
seplab = c('Sep04', 'Sep05', 'Sep06', 'Sep07', 'Sep08', 'Sep09', 'Sep10', 'Sep11')
seplabs = rep(seplab, 62)

rownames(sep_emas) = seplabs
colnames(sep_emas) = 'EMA'

#add in factor ID for the 8 levels 
sep04 = as.numeric(rownames(sep_emas) == 'Sep04')
sep05 = as.numeric(rownames(sep_emas) == 'Sep05')
sep06 = as.numeric(rownames(sep_emas) == 'Sep06')
sep07 = as.numeric(rownames(sep_emas) == 'Sep07')
sep08 = as.numeric(rownames(sep_emas) == 'Sep08')
sep09 = as.numeric(rownames(sep_emas) == 'Sep09')
sep10 = as.numeric(rownames(sep_emas) == 'Sep10')
sep11 = as.numeric(rownames(sep_emas) == 'Sep11')

sep_ema = as.numeric(sep_emas)
sep_ema_df = data.frame(sep_ema, sep04,sep05,sep06,sep07,sep08,sep09,sep10,sep11)

sep_ema_df$new[sep_ema_df$sep04==1]<-4
sep_ema_df$new[sep_ema_df$sep05==1]<-5
sep_ema_df$new[sep_ema_df$sep06==1]<-6
sep_ema_df$new[sep_ema_df$sep07==1]<-7
sep_ema_df$new[sep_ema_df$sep08==1]<-8
sep_ema_df$new[sep_ema_df$sep09==1]<-9
sep_ema_df$new[sep_ema_df$sep10==1]<-10
sep_ema_df$new[sep_ema_df$sep11==1]<-11


#September density plot  

sep.dens1 = ggplot(sep_ema_df, aes(x=sep_ema, group = new)) + 
  geom_density(aes(colour=factor(new)), size = 1, alpha = .60,kernel='epanechnikov') +
  xlab(NULL) + ylab(NULL) +
  scale_colour_manual(values = cbbPalette,
                      labels= c('`04','`05','`06','`07','`08','`09','`10','`11'),
                      name= '') +
                        scale_y_continuous(limits=c(0.0000, 0.0030), labels = percent) +
                        scale_x_continuous(limits=c(-2000,2000)) +
                        theme(axis.line = element_line(colour = "black"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.border = element_blank(),
                              panel.background = element_blank(),
                              legend.title=element_blank(),
                              legend.position = 'bottom', 
                              legend.text = element_text(color= 'black', size = 12, face = 'bold'),
                              legend.key = element_rect(fill = 'white', size = 0.5, linetype='solid'),
                              legend.key.size = unit(1.5, 'lines'))


print(sep.dens1)


setwd('c:/Users/peytonk/Dropbox/Pokies/Replication')
ggsave(plot=sep.dens1, filename=paste('sepdensV2', '.pdf', sep=''), scale=1,height=5, width=7)


#---------------------------------------October-------------------------------------------------------#
#create indicators for Oct04 ... Oct11, then attach to data frame w/ estimates
octlab = c('Oct04', 'Oct05', 'Oct06', 'Oct07', 'Oct08', 'Oct09', 'Oct10', 'Oct11')
octlabs = rep(octlab, 62)

rownames(oct_emas) = octlabs
colnames(oct_emas) = 'EMA'

#add in factor ID for the 8 levels 
oct04 = as.numeric(rownames(oct_emas) == 'Oct04')
oct05 = as.numeric(rownames(oct_emas) == 'Oct05')
oct06 = as.numeric(rownames(oct_emas) == 'Oct06')
oct07 = as.numeric(rownames(oct_emas) == 'Oct07')
oct08 = as.numeric(rownames(oct_emas) == 'Oct08')
oct09 = as.numeric(rownames(oct_emas) == 'Oct09')
oct10 = as.numeric(rownames(oct_emas) == 'Oct10')
oct11 = as.numeric(rownames(oct_emas) == 'Oct11')

oct_ema = as.numeric(oct_emas)
oct_ema_df = data.frame(oct_ema, oct04,oct05,oct06,oct07,oct08,oct09,oct10,oct11)

oct_ema_df$new[oct_ema_df$oct04==1]<-4
oct_ema_df$new[oct_ema_df$oct05==1]<-5
oct_ema_df$new[oct_ema_df$oct06==1]<-6
oct_ema_df$new[oct_ema_df$oct07==1]<-7
oct_ema_df$new[oct_ema_df$oct08==1]<-8
oct_ema_df$new[oct_ema_df$oct09==1]<-9
oct_ema_df$new[oct_ema_df$oct10==1]<-10
oct_ema_df$new[oct_ema_df$oct11==1]<-11


#October density plot  

oct.dens1 = ggplot(oct_ema_df, aes(x=oct_ema, group = new)) + 
  geom_density(aes(colour=factor(new)), size = 1, alpha = .60, kernel='epanechnikov') +
  xlab(NULL) + ylab(NULL) +
  scale_colour_manual(values = cbbPalette,
                      labels= c('`04','`05','`06','`07','`08','`09','`10','`11'),
                      name= '') +
                        scale_y_continuous(limits=c(0.0000, 0.0030), labels = percent) +
                        scale_x_continuous(limits=c(-2000,2000)) +
                        theme(axis.line = element_line(colour = "black"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.border = element_blank(),
                              panel.background = element_blank(),
                              legend.title=element_blank(),
                              legend.position = 'bottom', 
                              legend.text = element_text(color= 'black', size = 12, face = 'bold'),
                              legend.key = element_rect(fill = 'white', size = 0.5, linetype='solid'),
                              legend.key.size = unit(1.5, 'lines'))


print(oct.dens1)


setwd('c:/Users/peytonk/Dropbox/Pokies/Replication')
ggsave(plot=oct.dens1, filename=paste('octdensV2', '.pdf', sep=''), scale=1,height=5, width=7)




#---------------------------------------November-------------------------------------------------------#

#create indicators for Nov04 ... Nov11, then attach to data frame w/ estimates
novlab = c('Nov04', 'Nov05', 'Nov06', 'Nov07', 'Nov08', 'Nov09', 'Nov10', 'Nov11')
novlabs = rep(novlab, 62)

rownames(nov_emas) = novlabs
colnames(nov_emas) = 'EMA'

#add in factor ID for the 8 levels 
nov04 = as.numeric(rownames(nov_emas) == 'Nov04')
nov05 = as.numeric(rownames(nov_emas) == 'Nov05')
nov06 = as.numeric(rownames(nov_emas) == 'Nov06')
nov07 = as.numeric(rownames(nov_emas) == 'Nov07')
nov08 = as.numeric(rownames(nov_emas) == 'Nov08')
nov09 = as.numeric(rownames(nov_emas) == 'Nov09')
nov10 = as.numeric(rownames(nov_emas) == 'Nov10')
nov11 = as.numeric(rownames(nov_emas) == 'Nov11')

nov_ema = as.numeric(nov_emas)
nov_ema_df = data.frame(nov_ema, nov04,nov05,nov06,nov07,nov08,nov09,nov10,nov11)

nov_ema_df$new[nov_ema_df$nov04==1]<-4
nov_ema_df$new[nov_ema_df$nov05==1]<-5
nov_ema_df$new[nov_ema_df$nov06==1]<-6
nov_ema_df$new[nov_ema_df$nov07==1]<-7
nov_ema_df$new[nov_ema_df$nov08==1]<-8
nov_ema_df$new[nov_ema_df$nov09==1]<-9
nov_ema_df$new[nov_ema_df$nov10==1]<-10
nov_ema_df$new[nov_ema_df$nov11==1]<-11


#November density plot  
nov.dens1 = ggplot(nov_ema_df, aes(x=nov_ema, group = new)) + 
  geom_density(aes(colour=factor(new)), size = 1, alpha = .60,kernel='epanechnikov') +
  xlab(NULL) + ylab(NULL) +
  scale_colour_manual(values = cbbPalette,
                      labels= c('`04','`05','`06','`07','`08','`09','`10','`11'),
                      name= '') +
                        scale_y_continuous(limits=c(0.0000, 0.0030), labels = percent) +
                        scale_x_continuous(limits=c(-2000,2000)) +
                        theme(axis.line = element_line(colour = "black"),
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(),
                              panel.border = element_blank(),
                              panel.background = element_blank(),
                              legend.title=element_blank(),
                              legend.position = 'bottom', 
                              legend.text = element_text(color= 'black', size = 12, face = 'bold'),
                              legend.key = element_rect(fill = 'white', size = 0.5, linetype='solid'),
                              legend.key.size = unit(1.5, 'lines'))


print(nov.dens1)


setwd('c:/Users/peytonk/Dropbox/Pokies/Replication')
ggsave(plot=nov.dens1, filename=paste('novdensV2', '.pdf', sep=''), scale=1,height=5, width=7)

